{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "17ecfc5a",
      "metadata": {
        "id": "17ecfc5a"
      },
      "source": [
        "# Modeling Session 1\n",
        "This set of exercises are meant to cover many of the topics from Day 1, with callbacks from Opti 101, 201, 202 along the way. You can work along in this notebook, but if you really just want to see the solutions, those can also be found in the GitHub repo."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 51,
      "id": "hpOPsEkIbZ6E",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "hpOPsEkIbZ6E",
        "outputId": "a801df04-12ad-4e68-bb8d-ac11a4c15327"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Note: you may need to restart the kernel to use updated packages.\n"
          ]
        }
      ],
      "source": [
        "%pip -q install gurobipy==13.0.0"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 52,
      "id": "4449d5c7",
      "metadata": {
        "id": "4449d5c7"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import gurobipy as gp\n",
        "from gurobipy import GRB"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "bec77426",
      "metadata": {
        "id": "bec77426"
      },
      "source": [
        "## Building Trust in Optimization"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ebefb4a0",
      "metadata": {
        "id": "ebefb4a0"
      },
      "source": [
        "### Solution pools"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3880ef93",
      "metadata": {
        "id": "3880ef93"
      },
      "source": [
        "In this section, we'll explore solution pools to build trust in a discrete model by\n",
        "surfacing *multiple, near‑optimal* operating plans. We'll use the same storage + solar\n",
        "dispatch context from earlier (charging/discharging a battery, buying from the grid, meeting\n",
        "demand), but we’ll ask Gurobi to return a *set* of high‑quality alternatives.\n",
        "\n",
        "Goal: Collect the top‑k solutions within a small optimality gap and compare them on\n",
        "key features (objective, total grid energy, number of charge/discharge switches).\n",
        "\n",
        "Note that some of this (the problem, and the approaches) should seem familiar..."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "691fed82",
      "metadata": {},
      "source": [
        "#### 🔋 Energy Storage Optimization Model\n",
        "\n",
        "##### Sets\n",
        "$$\n",
        "\\begin{aligned}\n",
        "T = \\{0, 1, \\dots, 11\\} \\quad \\text{(time periods, e.g., hours)}\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "##### Parameters\n",
        "$$\n",
        "\\begin{aligned}\n",
        "p_t & : \\text{Electricity price in period } t \\; [€/MWh] \\\\\n",
        "d_t & : \\text{Energy demand in period } t \\; [\\text{MWh}] \\\\\n",
        "\\eta_{\\text{in}}, \\eta_{\\text{out}} & : \\text{Charging and discharging efficiencies} \\\\\n",
        "\\text{cap} & : \\text{Battery capacity } [\\text{MWh}] \\\\\n",
        "p_{\\max} & : \\text{Maximum charge/discharge power } [\\text{MW}] \\\\\n",
        "s_0 & : \\text{Initial state of charge } [\\text{MWh}]\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "##### Decision Variables\n",
        "$$\n",
        "\\begin{aligned}\n",
        "ch_t &\\ge 0 && \\text{Charging power (MW)} \\\\\n",
        "dis_t &\\ge 0 && \\text{Discharging power (MW)} \\\\\n",
        "grid_t &\\ge 0 && \\text{Power purchased from the grid (MW)} \\\\\n",
        "0 \\le soc_t &\\le \\text{cap} && \\text{State of charge (MWh)} \\\\\n",
        "y_{ch,t} &\\in \\{0,1\\} && \\text{1 if charging in period } t \\\\\n",
        "y_{dis,t} &\\in \\{0,1\\} && \\text{1 if discharging in period } t\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "##### Objective Function\n",
        "Minimize total grid cost and a small switching penalty:\n",
        "$$\n",
        "\\begin{aligned}\n",
        "\\min \\; Z = \\sum_{t \\in T} p_t \\cdot grid_t \\;+\\; 10^{-3} \\sum_{t \\in T} (y_{ch,t} + y_{dis,t})\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "\n",
        "##### Constraints\n",
        "\n",
        "1. Mutual Exclusivity\n",
        "$$\n",
        "\\begin{aligned}\n",
        "y_{ch,t} + y_{dis,t} \\le 1 \\quad \\forall t \\in T\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "2. Power Limits\n",
        "$$\n",
        "\\begin{aligned}\n",
        "ch_t &\\le p_{\\max} \\, y_{ch,t} \\\\\n",
        "dis_t &\\le p_{\\max} \\, y_{dis,t} \\quad \\forall t \\in T\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "3. Energy Balance\n",
        "$$\n",
        "\\begin{aligned}\n",
        "grid_t + \\eta_{\\text{out}} \\, dis_t \\ge d_t + ch_t \\quad \\forall t \\in T\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "4. State-of-Charge (SoC) Dynamics\n",
        "$$\n",
        "\\begin{aligned}\n",
        "soc_0 &= s_0 + \\eta_{\\text{in}} \\, ch_0 - \\frac{1}{\\eta_{\\text{out}}} \\, dis_0 \\\\\n",
        "soc_t &= soc_{t-1} + \\eta_{\\text{in}} \\, ch_t - \\frac{1}{\\eta_{\\text{out}}} \\, dis_t \\quad \\forall t > 0 \\\\\n",
        "&0 \\le soc_t \\le \\text{cap} \\quad \\forall t \\in T\n",
        "\\end{aligned}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "aaa99a3d",
      "metadata": {},
      "source": [
        "For reference, here is our [solution pool documentation](https://docs.gurobi.com/projects/optimizer/en/current/features/solutionpool.html)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 53,
      "id": "20c28f6a",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "20c28f6a",
        "outputId": "b233bad9-1fbc-46d2-ef83-69e6d99fa0f5"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter PoolSearchMode to value 2\n",
            "Set parameter PoolSolutions to value 500\n",
            "Set parameter PoolGap to value 0.02\n",
            "Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[rosetta2] - Darwin 24.6.0 24G325)\n",
            "\n",
            "CPU model: Apple M4\n",
            "Thread count: 10 physical cores, 10 logical processors, using up to 10 threads\n",
            "\n",
            "Non-default parameters:\n",
            "PoolSolutions  500\n",
            "PoolSearchMode  2\n",
            "PoolGap  0.02\n",
            "\n",
            "Optimize a model with 60 rows, 72 columns and 155 nonzeros (Min)\n",
            "Model fingerprint: 0xd2c985b1\n",
            "Model has 36 linear objective coefficients\n",
            "Variable types: 48 continuous, 24 integer (24 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [9e-01, 4e+00]\n",
            "  Objective range  [1e-03, 4e+01]\n",
            "  Bounds range     [1e+00, 1e+01]\n",
            "  RHS range        [1e+00, 7e+00]\n",
            "Found heuristic solution: objective 1770.9009999\n",
            "Presolve removed 1 rows and 1 columns\n",
            "Presolve time: 0.00s\n",
            "Presolved: 59 rows, 71 columns, 153 nonzeros\n",
            "Variable types: 47 continuous, 24 integer (24 binary)\n",
            "\n",
            "Root relaxation: objective 1.562452e+03, 48 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 1562.45206    0    2 1770.90100 1562.45206  11.8%     -    0s\n",
            "H    0     0                    1562.4530526 1562.45206  0.00%     -    0s\n",
            "     0     0          -    0      1562.45305 1562.45305  0.00%     -    0s\n",
            "\n",
            "Optimal solution found at node 0 - now completing solution pool...\n",
            "\n",
            "    Nodes    |    Current Node    |      Pool Obj. Bounds     |     Work\n",
            "             |                    |   Worst                   |\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0          -    0               - 1562.45305      -     -    0s\n",
            "     0     0          -    0               - 1562.45305      -     -    0s\n",
            "     0     0          -    0               - 1562.45305      -     -    0s\n",
            "     0     2          -    0               - 1562.45305      -     -    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Gomory: 2\n",
            "  MIR: 4\n",
            "  Flow cover: 2\n",
            "\n",
            "Explored 3875 nodes (2689 simplex iterations) in 0.04 seconds (0.01 work units)\n",
            "Thread count was 10 (of 10 available processors)\n",
            "\n",
            "Solution count 500: 1562.45 1562.45 1562.45 ... 1564.87\n",
            "No other solutions better than 1564.87\n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 1.562453052632e+03, best bound 1.562453052632e+03, gap 0.0000%\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>obj</th>\n",
              "      <th>grid_MWh</th>\n",
              "      <th>switch_count</th>\n",
              "      <th>rel_diff_to_best_pct</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>1562.453053</td>\n",
              "      <td>69.696316</td>\n",
              "      <td>7.0</td>\n",
              "      <td>0.0%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>1562.454053</td>\n",
              "      <td>69.696316</td>\n",
              "      <td>8.0</td>\n",
              "      <td>6.4e-05%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>1562.455053</td>\n",
              "      <td>69.696316</td>\n",
              "      <td>9.0</td>\n",
              "      <td>0.000128%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>1562.456053</td>\n",
              "      <td>69.696316</td>\n",
              "      <td>10.0</td>\n",
              "      <td>0.000192%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>1562.457053</td>\n",
              "      <td>69.696316</td>\n",
              "      <td>11.0</td>\n",
              "      <td>0.000256%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>1562.458053</td>\n",
              "      <td>69.696316</td>\n",
              "      <td>12.0</td>\n",
              "      <td>0.00032%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>1564.867842</td>\n",
              "      <td>69.396053</td>\n",
              "      <td>6.0</td>\n",
              "      <td>0.154551%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>7</th>\n",
              "      <td>1564.868842</td>\n",
              "      <td>69.396053</td>\n",
              "      <td>7.0</td>\n",
              "      <td>0.154615%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>8</th>\n",
              "      <td>1564.869842</td>\n",
              "      <td>69.396053</td>\n",
              "      <td>8.0</td>\n",
              "      <td>0.154679%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>9</th>\n",
              "      <td>1564.870842</td>\n",
              "      <td>69.396053</td>\n",
              "      <td>9.0</td>\n",
              "      <td>0.154743%</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>10</th>\n",
              "      <td>1564.871842</td>\n",
              "      <td>69.396053</td>\n",
              "      <td>10.0</td>\n",
              "      <td>0.154807%</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "            obj   grid_MWh  switch_count rel_diff_to_best_pct\n",
              "0   1562.453053  69.696316           7.0                 0.0%\n",
              "1   1562.454053  69.696316           8.0             6.4e-05%\n",
              "2   1562.455053  69.696316           9.0            0.000128%\n",
              "3   1562.456053  69.696316          10.0            0.000192%\n",
              "4   1562.457053  69.696316          11.0            0.000256%\n",
              "5   1562.458053  69.696316          12.0             0.00032%\n",
              "6   1564.867842  69.396053           6.0            0.154551%\n",
              "7   1564.868842  69.396053           7.0            0.154615%\n",
              "8   1564.869842  69.396053           8.0            0.154679%\n",
              "9   1564.870842  69.396053           9.0            0.154743%\n",
              "10  1564.871842  69.396053          10.0            0.154807%"
            ]
          },
          "execution_count": 53,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "import pandas as pd\n",
        "import gurobipy as gp\n",
        "from gurobipy import GRB\n",
        "\n",
        "# Tiny instance (kept small so we can enumerate several pool solutions quickly)\n",
        "T = range(12)                        # time periods\n",
        "price = [20,18,16,15,16,18,25,40,45,35,28,22]  # grid price €/MWh\n",
        "demand = [5,5,5,5,6,6,7,7,7,6,6,5]  # MWh demand each period\n",
        "\n",
        "eta_in, eta_out = 0.95, 0.95        # round-trip efficiency\n",
        "cap = 10.0                          # MWh battery capacity\n",
        "pmax = 4.0                          # MW charge/discharge limit\n",
        "s0 = 2.0                            # initial state of charge\n",
        "\n",
        "m = gp.Model(\"storage_pool\")\n",
        "\n",
        "# Decision variables\n",
        "ch = m.addVars(T, name=\"charge\", lb=0, ub=pmax)\n",
        "dis = m.addVars(T, name=\"discharge\", lb=0, ub=pmax)\n",
        "grid = m.addVars(T, name=\"grid\", lb=0)\n",
        "soc = m.addVars(T, name=\"soc\", lb=0, ub=cap)\n",
        "# on/off switching to discourage rapid toggling (makes model MIP)\n",
        "y_ch = m.addVars(T, vtype=GRB.BINARY, name=\"y_ch\")\n",
        "y_dis = m.addVars(T, vtype=GRB.BINARY, name=\"y_dis\")\n",
        "\n",
        "# Logic: can’t charge and discharge simultaneously\n",
        "m.addConstrs((y_ch[t] + y_dis[t] <= 1 for t in T), name=\"no_simul\")\n",
        "m.addConstrs((ch[t] <= pmax * y_ch[t] for t in T), name=\"charge_gate\")\n",
        "m.addConstrs((dis[t] <= pmax * y_dis[t] for t in T), name=\"discharge_gate\")\n",
        "\n",
        "# Energy balance vs demand\n",
        "m.addConstrs((grid[t] + dis[t]*eta_out >= demand[t] + ch[t] for t in T), name=\"balance\")\n",
        "\n",
        "# State of charge dynamics\n",
        "for t in T:\n",
        "    if t == 0:\n",
        "        m.addConstr(soc[t] == s0 + ch[t]*eta_in - dis[t]/eta_out, name=f\"soc_{t}\")\n",
        "    else:\n",
        "        m.addConstr(soc[t] == soc[t-1] + ch[t]*eta_in - dis[t]/eta_out, name=f\"soc_{t}\")\n",
        "\n",
        "# Mild cycling preference: penalize on/off changes to diversify plans\n",
        "switches = gp.quicksum((y_ch[t] + y_dis[t]) for t in T)\n",
        "\n",
        "# Cost: energy purchased from grid\n",
        "cost = gp.quicksum(price[t]*grid[t] for t in T)\n",
        "\n",
        "# Multi-criteria via a single objective (lexico done elsewhere; here we pool around cost)\n",
        "m.setObjective(cost + 1e-3 * switches, GRB.MINIMIZE)\n",
        "\n",
        "# Ask for a pool of near‑optimal solutions\n",
        "m.setParam(GRB.Param.PoolSearchMode, 2)  # systematic search for k‑best\n",
        "m.setParam(GRB.Param.PoolSolutions, 500)  # keep up to 500\n",
        "m.setParam(GRB.Param.PoolGap, 0.02)      # within 2% of best found\n",
        "m.optimize()\n",
        "\n",
        "# Collect solutions\n",
        "rows = []\n",
        "for i in range(m.SolCount):\n",
        "    m.setParam(GRB.Param.SolutionNumber, i)\n",
        "    obj = m.PoolObjVal\n",
        "    total_grid = sum(grid[t].Xn for t in T)\n",
        "    total_switch = sum(y_ch[t].Xn + y_dis[t].Xn for t in T)\n",
        "    rows.append({\"obj\": obj, \"grid_MWh\": total_grid, \"switch_count\": total_switch})\n",
        "\n",
        "sols = pd.DataFrame(rows).sort_values(\"obj\").reset_index(drop=True)\n",
        "\n",
        "best_obj = m.ObjVal\n",
        "sols[\"rel_diff_to_best_pct\"] = ((sols[\"obj\"] - best_obj) / abs(best_obj)* 100).round(6).astype(str) + \"%\"\n",
        "sols_unique = sols.drop_duplicates().reset_index(drop=True)\n",
        "sols_unique\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8ac4f028",
      "metadata": {},
      "source": [
        "<div style=\"\n",
        "  border-left: 6px solid #f59e0b;  /* amber/yellow */\n",
        "  padding:14px 16px;\n",
        "  border-radius:10px;\n",
        "  margin:12px 0;\n",
        "  font-size:16px; line-height:1.4;\">\n",
        "  <strong style=\"font-size:18px;\">💬 Explainability Moment</strong><br>\n",
        "  <em>Pause here.</em> How would you explain the methodology of solution pools and this group of solutions to a <b>non-technical stakeholder</b>? \n",
        "</div>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "184ff49f",
      "metadata": {
        "id": "184ff49f"
      },
      "source": [
        "### Multiple scenarios"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "16e56094",
      "metadata": {
        "id": "16e56094"
      },
      "source": [
        "In this section, we'll evaluate multiple scenarios by changing the *demand* and *price*\n",
        "profiles and re‑solving the same dispatch model. This mirrors the “what‑if” analysis you saw\n",
        "earlier: same structure, new data.\n",
        "\n",
        "Goal: Run three scenarios (Base, High‑Price Peak, Evening‑Heavy Demand) and compare total\n",
        "grid cost and final state of charge."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 54,
      "id": "f9b59570",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 143
        },
        "id": "f9b59570",
        "outputId": "d2dca684-87d6-4946-a00c-5c0e2ff8df1b"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>scenario</th>\n",
              "      <th>total_cost</th>\n",
              "      <th>final_soc</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>Base</td>\n",
              "      <td>1562.45</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>HighPricePeak</td>\n",
              "      <td>1648.87</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>EveningHeavyDemand</td>\n",
              "      <td>1779.05</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "             scenario  total_cost  final_soc\n",
              "0                Base     1562.45        0.0\n",
              "1       HighPricePeak     1648.87        0.0\n",
              "2  EveningHeavyDemand     1779.05        0.0"
            ]
          },
          "execution_count": 54,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "import pandas as pd\n",
        "import gurobipy as gp\n",
        "from gurobipy import GRB\n",
        "\n",
        "T = range(12)\n",
        "base_price = [20,18,16,15,16,18,25,40,45,35,28,22]\n",
        "base_demand = [5,5,5,5,6,6,7,7,7,6,6,5]\n",
        "\n",
        "scenarios = {\n",
        "    \"Base\": (base_price, base_demand),\n",
        "    \"HighPricePeak\": ([p*1.2 if 7<=t<=9 else p for t,p in enumerate(base_price)], base_demand),\n",
        "    \"EveningHeavyDemand\": (base_price, [d*1.2 if t>=7 else d for t,d in enumerate(base_demand)]),\n",
        "}\n",
        "\n",
        "eta_in, eta_out = 0.95, 0.95\n",
        "cap, pmax, s0 = 10.0, 4.0, 2.0\n",
        "\n",
        "def solve(price, demand):\n",
        "    m = gp.Model(\"storage_scen\")\n",
        "    ch = m.addVars(T, name=\"charge\", lb=0, ub=pmax)\n",
        "    dis = m.addVars(T, name=\"discharge\", lb=0, ub=pmax)\n",
        "    grid = m.addVars(T, name=\"grid\", lb=0)\n",
        "    soc = m.addVars(T, name=\"soc\", lb=0, ub=cap)\n",
        "\n",
        "    # LP version to keep focus on param changes\n",
        "    m.addConstrs((grid[t] + dis[t]*eta_out >= demand[t] + ch[t] for t in T), name=\"balance\")\n",
        "    for t in T:\n",
        "        if t == 0:\n",
        "            m.addConstr(soc[t] == s0 + ch[t]*eta_in - dis[t]/eta_out)\n",
        "        else:\n",
        "            m.addConstr(soc[t] == soc[t-1] + ch[t]*eta_in - dis[t]/eta_out)\n",
        "\n",
        "    cost = gp.quicksum(price[t]*grid[t] for t in T)\n",
        "    m.setObjective(cost, GRB.MINIMIZE)\n",
        "    m.setParam(\"OutputFlag\", 0)\n",
        "    m.optimize()\n",
        "\n",
        "    return float(cost.getValue()), float(soc[max(T)].X)\n",
        "\n",
        "rows = []\n",
        "for name, (p, d) in scenarios.items():\n",
        "    total_cost, final_soc = solve(p, d)\n",
        "    rows.append({\"scenario\": name, \"total_cost\": round(total_cost,2), \"final_soc\": round(final_soc,2)})\n",
        "\n",
        "pd.DataFrame(rows).sort_values(\"total_cost\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6ea8c166",
      "metadata": {},
      "source": [
        "<div style=\"\n",
        "  border-left: 6px solid #f59e0b;  /* amber/yellow */\n",
        "  padding:14px 16px;\n",
        "  border-radius:10px;\n",
        "  margin:12px 0;\n",
        "  font-size:16px; line-height:1.4;\">\n",
        "  <strong style=\"font-size:18px;\">💬 Explainability Moment</strong><br>\n",
        "  <em>Pause here.</em> How do you think running multiple scenarios can build trust in your optimization model? What other scenarios do you think could be helpful to run for this problem?\n",
        "</div>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3e07dc1e",
      "metadata": {},
      "source": [
        "<div style=\"\n",
        "  border-left: 6px solid #3b82f6;  /* blue */\n",
        "  padding:14px 16px;\n",
        "  border-radius:10px;\n",
        "  margin:12px 0;\n",
        "  font-size:16px; line-height:1.4;\">\n",
        "  <strong style=\"font-size:18px;\">📘 Optional Activity: Back to Hidden Gems</strong><br>\n",
        "  Using Gurobi's multi-scenario functionality to solve the same problem above. \n",
        "</div>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "408732ff",
      "metadata": {
        "id": "408732ff"
      },
      "source": [
        "### Pi values and sensitivity analysis"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "976977ed",
      "metadata": {
        "id": "976977ed"
      },
      "source": [
        "In this section, we'll inspect dual values (π) and sensitivity for the *LP* version\n",
        "of the dispatch model. The dual (shadow) price on a constraint tells you how much the\n",
        "objective would improve per unit relaxation of that constraint’s RHS - holding the basis fixed.\n",
        "\n",
        "Goal: Solve the LP, read π for the balance constraints, and identify which periods\n",
        "are most “binding” (largest |π|). Then, perturb demand slightly and see if the objective\n",
        "change agrees with π."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 55,
      "id": "571b6297",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 319
        },
        "id": "571b6297",
        "outputId": "3ce404cd-81e3-4969-cd96-faaf3054da2b"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Top periods by shadow price (higher pi => relaxing balance helps more):\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>t</th>\n",
              "      <th>pi</th>\n",
              "      <th>slack</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>8</th>\n",
              "      <td>8</td>\n",
              "      <td>45.0</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>7</th>\n",
              "      <td>7</td>\n",
              "      <td>40.0</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>9</th>\n",
              "      <td>9</td>\n",
              "      <td>35.0</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>10</th>\n",
              "      <td>10</td>\n",
              "      <td>28.0</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>6</td>\n",
              "      <td>25.0</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>11</th>\n",
              "      <td>11</td>\n",
              "      <td>22.0</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "     t    pi  slack\n",
              "8    8  45.0    0.0\n",
              "7    7  40.0    0.0\n",
              "9    9  35.0    0.0\n",
              "10  10  28.0    0.0\n",
              "6    6  25.0    0.0\n",
              "11  11  22.0    0.0"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>t*</th>\n",
              "      <th>pi</th>\n",
              "      <th>bump</th>\n",
              "      <th>approx_obj_change (pi*bump)</th>\n",
              "      <th>actual_obj_change</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>8</td>\n",
              "      <td>45.0</td>\n",
              "      <td>0.1</td>\n",
              "      <td>4.5</td>\n",
              "      <td>4.5</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   t*    pi  bump  approx_obj_change (pi*bump)  actual_obj_change\n",
              "0   8  45.0   0.1                          4.5                4.5"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "FYI,\n",
            "- Period 8 had the highest shadow price (pi = 45.0), meaning extra demand there is most expensive.\n",
            "- The bump of 0.1 MWh is small enough that the linear, pi-based prediction is valid.\n",
            "- Predicted vs actual change (4.50 vs 4.50) are nearly identical, so we are within the linear sensitivity range.\n",
            "- Slack = 0 confirms this constraint was binding.\n",
            "- In plain terms, adding 0.1 MWh at hour 8 costs about €4.5, exactly as the model’s dual value predicted.\n"
          ]
        }
      ],
      "source": [
        "import pandas as pd\n",
        "import gurobipy as gp\n",
        "from gurobipy import GRB\n",
        "\n",
        "T = range(12)\n",
        "price = [20,18,16,15,16,18,25,40,45,35,28,22]\n",
        "demand = [5,5,5,5,6,6,7,7,7,6,6,5]\n",
        "\n",
        "eta_in, eta_out = 0.95, 0.95\n",
        "cap, pmax, s0 = 10.0, 4.0, 2.0\n",
        "\n",
        "# Build LP (no binaries)\n",
        "m = gp.Model(\"storage_duals\")\n",
        "ch = m.addVars(T, name=\"charge\", lb=0, ub=pmax)\n",
        "dis = m.addVars(T, name=\"discharge\", lb=0, ub=pmax)\n",
        "grid = m.addVars(T, name=\"grid\", lb=0)\n",
        "soc = m.addVars(T, name=\"soc\", lb=0, ub=cap)\n",
        "\n",
        "balance = m.addConstrs((grid[t] + dis[t]*eta_out >= demand[t] + ch[t] for t in T), name=\"balance\")\n",
        "for t in T:\n",
        "    if t == 0:\n",
        "        m.addConstr(soc[t] == s0 + ch[t]*eta_in - dis[t]/eta_out, name=f\"soc_{t}\")\n",
        "    else:\n",
        "        m.addConstr(soc[t] == soc[t-1] + ch[t]*eta_in - dis[t]/eta_out, name=f\"soc_{t}\")\n",
        "\n",
        "cost = gp.quicksum(price[t]*grid[t] for t in T)\n",
        "m.setObjective(cost, GRB.MINIMIZE)\n",
        "m.setParam(\"OutputFlag\", 0)\n",
        "m.optimize()\n",
        "\n",
        "# Read duals (Pi) and slacks\n",
        "rows = []\n",
        "for t in T:\n",
        "    con = balance[t]\n",
        "    rows.append({\n",
        "        \"t\": t,\n",
        "        \"pi\": con.Pi,\n",
        "        \"slack\": con.Slack,\n",
        "    })\n",
        "\n",
        "pi_df = pd.DataFrame(rows).sort_values(\"pi\", ascending=False)\n",
        "print(\"Top periods by shadow price (higher pi => relaxing balance helps more):\")\n",
        "display(pi_df.head(6))\n",
        "\n",
        "# Sensitivity check: bump demand at period with largest pi by +0.1 and compare objective change\n",
        "t_star = int(pi_df.iloc[0][\"t\"])\n",
        "bump = 0.1\n",
        "\n",
        "def solve_with_demand(d):\n",
        "    m2 = gp.Model(\"storage_duals_copy\")\n",
        "    ch2 = m2.addVars(T, name=\"charge\", lb=0, ub=pmax)\n",
        "    dis2 = m2.addVars(T, name=\"discharge\", lb=0, ub=pmax)\n",
        "    grid2 = m2.addVars(T, name=\"grid\", lb=0)\n",
        "    soc2 = m2.addVars(T, name=\"soc\", lb=0, ub=cap)\n",
        "    bal2 = m2.addConstrs((grid2[t] + dis2[t]*eta_out >= d[t] + ch2[t] for t in T))\n",
        "    for t in T:\n",
        "        if t == 0:\n",
        "            m2.addConstr(soc2[t] == s0 + ch2[t]*eta_in - dis2[t]/eta_out)\n",
        "        else:\n",
        "            m2.addConstr(soc2[t] == soc2[t-1] + ch2[t]*eta_in - dis2[t]/eta_out)\n",
        "    cost2 = gp.quicksum(price[t]*grid2[t] for t in T)\n",
        "    m2.setObjective(cost2, GRB.MINIMIZE)\n",
        "    m2.setParam(\"OutputFlag\", 0)\n",
        "    m2.optimize()\n",
        "    return float(cost2.getValue())\n",
        "\n",
        "base_obj = float(cost.getValue())\n",
        "d2 = demand.copy()\n",
        "d2[t_star] += bump\n",
        "new_obj = solve_with_demand(d2)\n",
        "\n",
        "approx_change = pi_df.iloc[0][\"pi\"] * bump\n",
        "actual_change = new_obj - base_obj\n",
        "\n",
        "summary = pd.DataFrame({\n",
        "    \"t*\":[t_star],\n",
        "    \"pi\":[pi_df.iloc[0][\"pi\"]],\n",
        "    \"bump\":[bump],\n",
        "    \"approx_obj_change (pi*bump)\":[approx_change],\n",
        "    \"actual_obj_change\":[actual_change],\n",
        "})\n",
        "display(summary)\n",
        " \n",
        "# FYI explanation printout, commas only\n",
        "print(\n",
        "    f\"\\nFYI,\\n\"\n",
        "    f\"- Period {t_star} had the highest shadow price (pi = {pi_df.iloc[0]['pi']:.1f}), meaning extra demand there is most expensive.\\n\"\n",
        "    f\"- The bump of {bump} MWh is small enough that the linear, pi-based prediction is valid.\\n\"\n",
        "    f\"- Predicted vs actual change ({approx_change:.2f} vs {actual_change:.2f}) are nearly identical, so we are within the linear sensitivity range.\\n\"\n",
        "    f\"- Slack = 0 confirms this constraint was binding.\\n\"\n",
        "    f\"- In plain terms, adding 0.1 MWh at hour {t_star} costs about €{approx_change:.1f}, exactly as the model’s dual value predicted.\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2f5c6372",
      "metadata": {},
      "source": [
        "### Cautionary Tale: When Duals Lie to You\n",
        " \n",
        "Shadow prices (π’s) are powerful, but they only tell the truth *locally*, right around the current optimal solution, assuming nothing else changes.\n",
        " \n",
        "Imagine this:\n",
        " \n",
        "> The model says π = €45 at hour 8,  \n",
        "> you say “Adding 1 MWh will cost €45!”,  \n",
        "> you rerun the model, and the cost jumps by €70.\n",
        " \n",
        "**What happened?**\n",
        " \n",
        "- You crossed a kink in the feasible region, a new constraint became binding, so the old π no longer applies.  \n",
        "- You took too large a step, duals are like instantaneous slopes, not good for big jumps.  \n",
        "- You have integer logic (on/off decisions), duals only exist for the continuous relaxation, not for the real MIP behavior.  \n",
        "- Or your scaling or units are inconsistent, a π in kWh means something different than in MWh."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "28e3aa1b",
      "metadata": {},
      "source": [
        "<div style=\"\n",
        "  border-left: 6px solid #f59e0b;  /* amber/yellow */\n",
        "  padding:14px 16px;\n",
        "  border-radius:10px;\n",
        "  margin:12px 0;\n",
        "  font-size:16px; line-height:1.4;\">\n",
        "  <strong style=\"font-size:18px;\">💬 Moral of the story</strong><br>\n",
        "Duals are great for intuition and small “what-if” nudges, but risky for big decisions. However, they are like weather forecasts, useful for local conditions, unreliable once you cross into another climate.\n",
        "</div>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "fa98af39",
      "metadata": {
        "id": "fa98af39"
      },
      "source": [
        "## Optimization's Impact in Industry: Finance"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9a5b2948",
      "metadata": {
        "id": "9a5b2948"
      },
      "source": [
        "### Add a cardinality constraint to the basic Markowitz Portfolio Optimization model\n",
        "Now we'll apply some of the concepts of *discrete decisions* from the session on Optimization in Finance. A great place to start is by testing out the OptiMod for [Mean-Variance Portfolio Optimization](https://gurobi-optimods.readthedocs.io/en/stable/index.html). Then, we'll build our own version of the same model to compare. Finally, we'll change and add to that model, similar to what's shown in the [Gurobi Finance](https://gurobi-finance.readthedocs.io/en/latest/index.html) documentation, to reflect different scenarios."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 56,
      "id": "c369fde2",
      "metadata": {
        "id": "c369fde2"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Note: you may need to restart the kernel to use updated packages.\n"
          ]
        }
      ],
      "source": [
        "# Install the optimods package\n",
        "%pip -q install gurobi-optimods\n",
        "import numpy as np\n",
        "\n",
        "from gurobi_optimods.datasets import load_portfolio\n",
        "from gurobi_optimods.portfolio import MeanVariancePortfolio"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "eb31f8cb",
      "metadata": {
        "id": "eb31f8cb"
      },
      "source": [
        "All of the individual OptiMods come with data to be used right out of the box for each pre-built model. Let's load that and briefly see what form the data is in. If you're using an OptiMod, pay attention to the format of the data to make sure any new data is in the same form."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 57,
      "id": "07d41a5f",
      "metadata": {},
      "outputs": [],
      "source": [
        "np.set_printoptions(legacy='1.21')\n",
        "\n",
        "data = load_portfolio()\n",
        "cov_matrix = data.cov()\n",
        "mu = data.mean()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 58,
      "id": "2806a2aa",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "AA    0.001602\n",
              "BB    0.014219\n",
              "CC   -0.013629\n",
              "DD    0.007160\n",
              "EE    0.001796\n",
              "FF   -0.004670\n",
              "GG    0.002771\n",
              "HH    0.005117\n",
              "II   -0.000745\n",
              "JJ    0.004155\n",
              "dtype: float64"
            ]
          },
          "execution_count": 58,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Show the expected returns for our potential 'assets'\n",
        "mu"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 59,
      "id": "54a6e57a",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>AA</th>\n",
              "      <th>BB</th>\n",
              "      <th>CC</th>\n",
              "      <th>DD</th>\n",
              "      <th>EE</th>\n",
              "      <th>FF</th>\n",
              "      <th>GG</th>\n",
              "      <th>HH</th>\n",
              "      <th>II</th>\n",
              "      <th>JJ</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>AA</th>\n",
              "      <td>0.000256</td>\n",
              "      <td>0.000199</td>\n",
              "      <td>0.000077</td>\n",
              "      <td>0.000069</td>\n",
              "      <td>0.000145</td>\n",
              "      <td>0.000152</td>\n",
              "      <td>0.000160</td>\n",
              "      <td>0.000082</td>\n",
              "      <td>0.000086</td>\n",
              "      <td>0.000242</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>BB</th>\n",
              "      <td>0.000199</td>\n",
              "      <td>0.002965</td>\n",
              "      <td>0.001186</td>\n",
              "      <td>0.000650</td>\n",
              "      <td>0.001006</td>\n",
              "      <td>0.000736</td>\n",
              "      <td>0.000586</td>\n",
              "      <td>0.000154</td>\n",
              "      <td>0.000309</td>\n",
              "      <td>0.001625</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>CC</th>\n",
              "      <td>0.000077</td>\n",
              "      <td>0.001186</td>\n",
              "      <td>0.002303</td>\n",
              "      <td>0.000310</td>\n",
              "      <td>0.000776</td>\n",
              "      <td>0.000767</td>\n",
              "      <td>0.000405</td>\n",
              "      <td>0.000156</td>\n",
              "      <td>0.000132</td>\n",
              "      <td>0.000749</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>DD</th>\n",
              "      <td>0.000069</td>\n",
              "      <td>0.000650</td>\n",
              "      <td>0.000310</td>\n",
              "      <td>0.000591</td>\n",
              "      <td>0.000381</td>\n",
              "      <td>0.000255</td>\n",
              "      <td>0.000195</td>\n",
              "      <td>0.000075</td>\n",
              "      <td>0.000213</td>\n",
              "      <td>0.000445</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>EE</th>\n",
              "      <td>0.000145</td>\n",
              "      <td>0.001006</td>\n",
              "      <td>0.000776</td>\n",
              "      <td>0.000381</td>\n",
              "      <td>0.001424</td>\n",
              "      <td>0.000507</td>\n",
              "      <td>0.000331</td>\n",
              "      <td>-0.000007</td>\n",
              "      <td>0.000282</td>\n",
              "      <td>0.000753</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>FF</th>\n",
              "      <td>0.000152</td>\n",
              "      <td>0.000736</td>\n",
              "      <td>0.000767</td>\n",
              "      <td>0.000255</td>\n",
              "      <td>0.000507</td>\n",
              "      <td>0.000903</td>\n",
              "      <td>0.000270</td>\n",
              "      <td>0.000046</td>\n",
              "      <td>0.000112</td>\n",
              "      <td>0.000584</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>GG</th>\n",
              "      <td>0.000160</td>\n",
              "      <td>0.000586</td>\n",
              "      <td>0.000405</td>\n",
              "      <td>0.000195</td>\n",
              "      <td>0.000331</td>\n",
              "      <td>0.000270</td>\n",
              "      <td>0.000528</td>\n",
              "      <td>0.000101</td>\n",
              "      <td>0.000181</td>\n",
              "      <td>0.000368</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>HH</th>\n",
              "      <td>0.000082</td>\n",
              "      <td>0.000154</td>\n",
              "      <td>0.000156</td>\n",
              "      <td>0.000075</td>\n",
              "      <td>-0.000007</td>\n",
              "      <td>0.000046</td>\n",
              "      <td>0.000101</td>\n",
              "      <td>0.000581</td>\n",
              "      <td>-0.000032</td>\n",
              "      <td>-0.000068</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>II</th>\n",
              "      <td>0.000086</td>\n",
              "      <td>0.000309</td>\n",
              "      <td>0.000132</td>\n",
              "      <td>0.000213</td>\n",
              "      <td>0.000282</td>\n",
              "      <td>0.000112</td>\n",
              "      <td>0.000181</td>\n",
              "      <td>-0.000032</td>\n",
              "      <td>0.000659</td>\n",
              "      <td>0.000335</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>JJ</th>\n",
              "      <td>0.000242</td>\n",
              "      <td>0.001625</td>\n",
              "      <td>0.000749</td>\n",
              "      <td>0.000445</td>\n",
              "      <td>0.000753</td>\n",
              "      <td>0.000584</td>\n",
              "      <td>0.000368</td>\n",
              "      <td>-0.000068</td>\n",
              "      <td>0.000335</td>\n",
              "      <td>0.001644</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "          AA        BB        CC        DD        EE        FF        GG  \\\n",
              "AA  0.000256  0.000199  0.000077  0.000069  0.000145  0.000152  0.000160   \n",
              "BB  0.000199  0.002965  0.001186  0.000650  0.001006  0.000736  0.000586   \n",
              "CC  0.000077  0.001186  0.002303  0.000310  0.000776  0.000767  0.000405   \n",
              "DD  0.000069  0.000650  0.000310  0.000591  0.000381  0.000255  0.000195   \n",
              "EE  0.000145  0.001006  0.000776  0.000381  0.001424  0.000507  0.000331   \n",
              "FF  0.000152  0.000736  0.000767  0.000255  0.000507  0.000903  0.000270   \n",
              "GG  0.000160  0.000586  0.000405  0.000195  0.000331  0.000270  0.000528   \n",
              "HH  0.000082  0.000154  0.000156  0.000075 -0.000007  0.000046  0.000101   \n",
              "II  0.000086  0.000309  0.000132  0.000213  0.000282  0.000112  0.000181   \n",
              "JJ  0.000242  0.001625  0.000749  0.000445  0.000753  0.000584  0.000368   \n",
              "\n",
              "          HH        II        JJ  \n",
              "AA  0.000082  0.000086  0.000242  \n",
              "BB  0.000154  0.000309  0.001625  \n",
              "CC  0.000156  0.000132  0.000749  \n",
              "DD  0.000075  0.000213  0.000445  \n",
              "EE -0.000007  0.000282  0.000753  \n",
              "FF  0.000046  0.000112  0.000584  \n",
              "GG  0.000101  0.000181  0.000368  \n",
              "HH  0.000581 -0.000032 -0.000068  \n",
              "II -0.000032  0.000659  0.000335  \n",
              "JJ -0.000068  0.000335  0.001644  "
            ]
          },
          "execution_count": 59,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Show the covariance matrix\n",
        "cov_matrix"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1948cdbc",
      "metadata": {},
      "source": [
        "The *cardinality constraint* we will implement first will just limit the number of positions we can take (i.e. the number of assets to invest in), which we'll initially set to 3. We also see the risk-aversion parameter, `gamma`($\\gamma$). The number used for these examples is not specifically set, so feel free to change it. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 60,
      "id": "a73d36ee",
      "metadata": {
        "id": "a73d36ee",
        "outputId": "8f407a70-6536-4a00-93d6-d7dd143cf058"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The portfolio's expected return is 0.00723\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "AA    0.000000\n",
              "BB    0.126590\n",
              "CC    0.000000\n",
              "DD    0.470867\n",
              "EE    0.000000\n",
              "FF    0.000000\n",
              "GG    0.000000\n",
              "HH    0.402543\n",
              "II    0.000000\n",
              "JJ    0.000000\n",
              "dtype: float64"
            ]
          },
          "execution_count": 60,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# We can define this risk-aversion coefficient\n",
        "gamma = 20\n",
        "# Max positions K\n",
        "K = 3\n",
        "\n",
        "# Fill in the above info for expected returns and the covariance, respectively\n",
        "mvp = MeanVariancePortfolio(mu, cov_matrix)\n",
        "\n",
        "# Put in gamma, max positions. If you want to see the logs, set verbose = True\n",
        "mp = mvp.efficient_portfolio(gamma, max_positions=K, verbose = False)\n",
        "\n",
        "print(f\"The portfolio's expected return is {round(mp.ret,5)}\")\n",
        "\n",
        "# The portfolio's allocation:\n",
        "mp.x"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "18716d24",
      "metadata": {},
      "source": [
        "Now it's time to build the same model ourselves. Some parts of the code will use Gurobi's Python matrix API. Here are some [basic examples](https://support.gurobi.com/hc/en-us/articles/17278438215313-Tutorial-Getting-Started-with-the-Gurobi-Python-API).\n",
        "\n",
        "As a reminder, this is the problem we want to model.\n",
        "\n",
        "Let $S$ be the set of assets available, $\\mu$ is the expected return, and $\\Sigma$ is the covariance matrix. $K$ is the max number of assets to invest in. \n",
        "$$\n",
        "\\begin{aligned}\n",
        "\\max_{x \\in \\mathbb{R}^n}\\quad \n",
        "& \\mu^\\top x \\;-\\; \\tfrac{\\gamma}{2}\\, x^\\top \\Sigma\\, x \\\\\n",
        "\\text{s.t.}\\quad \n",
        "& \\sum_{i \\in S} x_i = 1,\\\\\n",
        "& x_i \\le b_i, \\quad \\forall i \\in S \\\\\n",
        "& x_i \\ge l*b_i, \\quad \\forall i \\in S \\\\\n",
        "& \\sum_{i \\in S} b_i \\le K,\\\\\n",
        "& 0 \\le x_i \\le 1, \\quad \\forall i \\in S\n",
        "\\end{aligned}\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 61,
      "id": "8fa4cf77",
      "metadata": {
        "id": "8fa4cf77",
        "outputId": "72d5bf34-fd8d-4c79-dde9-a37406d1c53c"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[rosetta2] - Darwin 24.6.0 24G325)\n",
            "\n",
            "CPU model: Apple M4\n",
            "Thread count: 10 physical cores, 10 logical processors, using up to 10 threads\n",
            "\n",
            "Optimize a model with 22 rows, 20 columns and 60 nonzeros (Max)\n",
            "Model fingerprint: 0x1d83aebc\n",
            "Model has 10 linear objective coefficients\n",
            "Model has 55 quadratic objective terms\n",
            "Variable types: 10 continuous, 10 integer (10 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e-02, 1e+00]\n",
            "  Objective range  [7e-04, 1e-02]\n",
            "  QObjective range [3e-04, 7e-02]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [1e+00, 3e+00]\n",
            "Presolve time: 0.00s\n",
            "Presolved: 32 rows, 25 columns, 85 nonzeros\n",
            "Presolved model has 60 quadratic objective terms\n",
            "Variable types: 15 continuous, 10 integer (10 binary)\n",
            "Found heuristic solution: objective -0.0122823\n",
            "Found heuristic solution: objective -0.0009564\n",
            "Found heuristic solution: objective 0.0000923\n",
            "\n",
            "Root relaxation: objective 3.301807e-03, 39 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0    0.00330    0    2    0.00009    0.00330  3477%     -    0s\n",
            "H    0     0                       0.0029496    0.00330  11.9%     -    0s\n",
            "H    0     0                       0.0032908    0.00330  0.33%     -    0s\n",
            "\n",
            "Cutting planes:\n",
            "  RLT: 1\n",
            "\n",
            "Explored 1 nodes (39 simplex iterations) in 0.02 seconds (0.00 work units)\n",
            "Thread count was 10 (of 10 available processors)\n",
            "\n",
            "Solution count 5: 0.00329085 0.00294963 9.23185e-05 ... -0.0122823\n",
            "No other solutions better than 0.00329085\n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 3.290846152959e-03, best bound 3.290846152959e-03, gap 0.0000%\n"
          ]
        }
      ],
      "source": [
        "#### Now model it yourself\n",
        "gamma = 20\n",
        "l = 0.01\n",
        "# Create an empty optimization model\n",
        "m = gp.Model()\n",
        "\n",
        "# Add variables: x[i] denotes the proportion invested in stock i\n",
        "# 0 <= x[i] <= 1\n",
        "x = m.addMVar(len(mu), lb=0, ub=1, name=\"x\")\n",
        "b = m.addMVar(len(mu), vtype=gp.GRB.BINARY, name=\"b\")\n",
        "\n",
        "# Budget constraint: all investments sum up to 1\n",
        "m.addConstr(x.sum() == 1, name=\"Budget_Constraint\")\n",
        "\n",
        "# Link the continuous x to the binary b\n",
        "# If x > 0 then b = 1 (also remember the bounds of x)\n",
        "m.addConstr(x <= b, name=\"Indicator\")\n",
        "\n",
        "# Minimal position; see formula (2) above\n",
        "m.addConstr(x >= l * b, name=\"Minimal_Position\")\n",
        "# Cardinality constraint: at most K positions\n",
        "cardinality_constr = m.addConstr(b.sum() <= 3, \"Cardinality\")\n",
        "\n",
        "# Define objective function: Maximize expected utility\n",
        "m.setObjective(\n",
        "    mu.to_numpy() @ x - gamma / 2 * (x @ cov_matrix.to_numpy() @ x), gp.GRB.MAXIMIZE\n",
        ")\n",
        "m.optimize()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 62,
      "id": "14671723",
      "metadata": {
        "id": "14671723",
        "outputId": "afc23a8c-59f0-4d48-a1ba-1f35e2dcd9b1"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>OptiMod</th>\n",
              "      <th>Our Model</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>AA</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>BB</th>\n",
              "      <td>0.126590</td>\n",
              "      <td>0.126590</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>CC</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>DD</th>\n",
              "      <td>0.470867</td>\n",
              "      <td>0.470867</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>EE</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>FF</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>GG</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>HH</th>\n",
              "      <td>0.402543</td>\n",
              "      <td>0.402543</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>II</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>JJ</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "     OptiMod  Our Model\n",
              "AA  0.000000   0.000000\n",
              "BB  0.126590   0.126590\n",
              "CC  0.000000   0.000000\n",
              "DD  0.470867   0.470867\n",
              "EE  0.000000   0.000000\n",
              "FF  0.000000   0.000000\n",
              "GG  0.000000   0.000000\n",
              "HH  0.402543   0.402543\n",
              "II  0.000000   0.000000\n",
              "JJ  0.000000   0.000000"
            ]
          },
          "execution_count": 62,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "s2 = pd.Series(x.X, index= mp.x.index)\n",
        "s1 = pd.Series(mp.x)\n",
        "df = pd.concat([s2,s1], axis=1, ignore_index=True)\n",
        "df.columns = ['OptiMod','Our Model']\n",
        "df"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2b6672ed",
      "metadata": {},
      "source": [
        "### Move risk to a constraint\n",
        "In the models above, the objective represented a utility function that balanced risk and returns. Let's make a couple of changes to the model. \n",
        "\n",
        "- Model the risk as constraint, where the objective is now just to maximize return.\n",
        "- Incorporate fixed and variable transaction costs. \n",
        "- Remove the cardinality constraints from before \n",
        "\n",
        "Write the formulation where:\n",
        "$$\n",
        "\\begin{aligned}\n",
        "&\\sigma_0 \\space\\text{maximal variance for the portfolio} \\\\\n",
        "&c  \\space\\text{is the fixed transaction costs for any asset, relative to total investment value} \\\\\n",
        "&l  \\space\\text{is the lower bound on position size} \\\\\n",
        "&f_i \\space\\text{is the variable transaction fee for each asset relative to total investment value} \n",
        "\\end{aligned}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "90a23f5f",
      "metadata": {},
      "source": [
        "$$\n",
        "\\begin{aligned}\n",
        "\\max_{x \\in \\mathbb{R}^n}\\quad &\\mu^\\top x\\\\\n",
        "\\text{s.t.}\\quad \n",
        "&x^\\top \\Sigma\\, x \\le \\sigma_0\\\\\n",
        "& \\sum_{i \\in S} x_i + c\\sum_{i \\in S} b_i + \\sum_{i \\in S} f_ix_i= 1\\\\\n",
        "& x_i \\le b_i, \\quad \\forall i \\in S \\\\\n",
        "& x_i \\ge l*b_i, \\quad \\forall i \\in S \\\\\n",
        "& 0 \\le x_i \\le 1 \\quad \\forall i \\in S \\\\\n",
        "\n",
        "& b_i \\in \\{0,1\\} \\quad \\forall i \\in S\n",
        "\\end{aligned}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0cf8e74a",
      "metadata": {},
      "source": [
        "Use the data below to write code for the above model."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 63,
      "id": "ee00f964",
      "metadata": {},
      "outputs": [],
      "source": [
        "data = load_portfolio()\n",
        "cov_matrix = data.cov()\n",
        "mu = data.mean()\n",
        "\n",
        "# Values for the model parameters:\n",
        "V = 0.0025  # Maximal admissible variance (sigma^2)\n",
        "l = 0.001  # Minimal position size\n",
        "c = 0.0001  # Fixed transaction costs\n",
        "f = 0.001 * np.ones(mu.shape)  # Variable transaction fees"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 64,
      "id": "0bd48744",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[rosetta2] - Darwin 24.6.0 24G325)\n",
            "\n",
            "CPU model: Apple M4\n",
            "Thread count: 10 physical cores, 10 logical processors, using up to 10 threads\n",
            "\n",
            "Optimize a model with 21 rows, 20 columns and 60 nonzeros (Max)\n",
            "Model fingerprint: 0x3d88eddd\n",
            "Model has 10 linear objective coefficients\n",
            "Model has 1 quadratic constraint\n",
            "Variable types: 10 continuous, 10 integer (10 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e-04, 1e+00]\n",
            "  QMatrix range    [1e-05, 3e-03]\n",
            "  Objective range  [7e-04, 1e-02]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [1e+00, 1e+00]\n",
            "  QRHS range       [3e-03, 3e-03]\n",
            "Presolve time: 0.00s\n",
            "Presolved: 21 rows, 20 columns, 60 nonzeros\n",
            "Presolved model has 1 quadratic constraint(s)\n",
            "Variable types: 10 continuous, 10 integer (10 binary)\n",
            "\n",
            "Root relaxation: objective 1.390322e-02, 3 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0    0.01390    0    1          -    0.01390      -     -    0s\n",
            "H    0     0                       0.0122953    0.01390  13.1%     -    0s\n",
            "H    0     0                       0.0134646    0.01390  3.26%     -    0s\n",
            "H    0     0                       0.0134646    0.01390  3.26%     -    0s\n",
            "     0     0          -    0         0.01346    0.01347  0.00%     -    0s\n",
            "\n",
            "Explored 1 nodes (5 simplex iterations) in 0.01 seconds (0.00 work units)\n",
            "Thread count was 10 (of 10 available processors)\n",
            "\n",
            "Solution count 2: 0.0134646 0.0122953 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 1.346460588027e-02, best bound 1.346505105619e-02, gap 0.0033%\n"
          ]
        }
      ],
      "source": [
        "# Create an empty optimization model\n",
        "m = gp.Model()\n",
        "\n",
        "# Add variables: x[i] denotes the proportion invested in stock i\n",
        "x = m.addMVar(len(mu), lb=0, ub=1, name=\"x\")\n",
        "# Add variables: b[i]=1 if stock i is held, and b[i]=0 otherwise\n",
        "b = m.addMVar(len(mu), vtype=gp.GRB.BINARY, name=\"b\")\n",
        "# Upper bound on variance\n",
        "m.addConstr(x @ cov_matrix.to_numpy() @ x <= V, name=\"Variance\")\n",
        "# Force x to 0 if not traded; see formula (1) above\n",
        "m.addConstr(x <= b, name=\"Indicator\")\n",
        "# Minimal position; see formula (2) above\n",
        "m.addConstr(x >= l * b, name=\"Minimal_Position\")\n",
        "# Budget constraint: all investments, costs, and fees sum up to 1\n",
        "budget_constr = m.addConstr(\n",
        "    x.sum() + c*b.sum() + f @ x == 1, name=\"Budget_Constraint\"\n",
        ")\n",
        "# Define objective: Maximize expected return\n",
        "m.setObjective(mu.to_numpy() @ x, gp.GRB.MAXIMIZE)\n",
        "\n",
        "m.optimize()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 65,
      "id": "bb919952",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0    0.000000\n",
              "1    0.894386\n",
              "2    0.000000\n",
              "3    0.104415\n",
              "4    0.000000\n",
              "5    0.000000\n",
              "6    0.000000\n",
              "7    0.000000\n",
              "8    0.000000\n",
              "9    0.000000\n",
              "dtype: float64"
            ]
          },
          "execution_count": 65,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pd.Series(x.X)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1dacbf1d",
      "metadata": {
        "id": "1dacbf1d"
      },
      "source": [
        "## Interpreting Log Files"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d97e4910",
      "metadata": {},
      "source": [
        "### MIPFocus Parameter Analysis\n",
        "\n",
        "#### Introduction\n",
        "\n",
        "This section looks at the impact of the **MIPFocus** parameter on solving the `markshare_4_0.mps` model. You can find many models to tinker with in the [MIPLIB - Mixed-Integer Programming Library](https://miplib.zib.de/index.html)\n",
        "\n",
        "##### What is MIPFocus?\n",
        "\n",
        "The MIPFocus parameter allows you to modify the high-level solution strategy in Gurobi's MIP solver:\n",
        "\n",
        "- **MIPFocus = 0** (Default): Balances finding new feasible solutions and proving optimality\n",
        "- **MIPFocus = 1**: Focuses more on finding good quality feasible solutions quickly\n",
        "- **MIPFocus = 2**: Focuses more attention on proving optimality (assumes finding solutions is easy)\n",
        "- **MIPFocus = 3**: Focuses on improving the best objective bound\n",
        "\n",
        "In this analysis, we'll compare MIPFocus values 0, 1, and 2 on the glass4 model.\n",
        "\n",
        "##### Metrics to Analyze\n",
        "\n",
        "We will compare:\n",
        "- Solution time\n",
        "- Number of nodes explored\n",
        "- Number of simplex iterations\n",
        "- Final optimality gap\n",
        "- Objective value\n",
        "- Solution status"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 66,
      "id": "c441be31",
      "metadata": {},
      "outputs": [],
      "source": [
        "import gurobipy as gp\n",
        "from gurobipy import GRB\n",
        "import pandas as pd\n",
        "import time"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "98c5054a",
      "metadata": {},
      "source": [
        "#### Helper Function to Solve and Collect Statistics\n",
        "\n",
        "We'll create a helper function that:\n",
        "1. Loads the model\n",
        "2. Sets the MIPFocus parameter\n",
        "3. Solves the model\n",
        "4. Collects performance metrics"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9f1f184f",
      "metadata": {},
      "source": [
        "<div style=\"\n",
        "  border-left: 6px solid #ef4444;  /* red */\n",
        "  padding:14px 16px;\n",
        "  border-radius:10px;\n",
        "  margin:12px 0;\n",
        "  font-size:16px; line-height:1.4;\">\n",
        "  <strong style=\"font-size:18px; color:#b91c1c;\">❗ IMPORTANT</strong><br>\n",
        "  Download the mps file and unzip/extract it. <br>\n",
        "  You can find the file here:\n",
        "  <a href=\"https://miplib.zib.de/instance_details_markshare_4_0.html\" target=\"_blank\" rel=\"noopener noreferrer\">\n",
        "    MIPLIB — markshare_4_0 instance page\n",
        "  </a>.\n",
        "  <br><br>\n",
        "  Once you do this, if you downloaded this notebook and are running this locally make sure the file is in the same folder. <br>\n",
        "  If you are in Colab, click on the folder on the left part of your screen. Then drag the mps file there.\n",
        "</div>\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d4f1a862",
      "metadata": {},
      "outputs": [],
      "source": [
        "mps_file = 'markshare_4_0.mps' "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "641b0dab",
      "metadata": {},
      "outputs": [],
      "source": [
        "def solve_with_mipfocus(model_file = mps_file, mipfocus_value = 0, output_flag = 1, time_limit=60):\n",
        "    \"\"\"\n",
        "    Solve a model with a specific MIPFocus value and return performance metrics.\n",
        "    \n",
        "    Parameters:\n",
        "    -----------\n",
        "    model_file : str\n",
        "        Path to the .mps file\n",
        "    mipfocus_value : int\n",
        "        MIPFocus parameter value (0, 1, 2, or 3)\n",
        "    time_limit : float\n",
        "        Time limit in seconds\n",
        "    \n",
        "    Returns:\n",
        "    --------\n",
        "    dict : Dictionary containing performance metrics\n",
        "    \"\"\"\n",
        "    # Create environment and model\n",
        "    # This will make the name of a '*.log' file\n",
        "    log_file = f'mipfocus{mipfocus_value}.log'\n",
        "    env = gp.Env(params={\"LogFile\": log_file, \n",
        "                         \"OutputFlag\": output_flag, \n",
        "                         \"LogToConsole\": 0})  \n",
        "    model = gp.read(model_file, env=env)\n",
        "    \n",
        "    # Set parameters\n",
        "    model.setParam('MIPFocus', mipfocus_value)\n",
        "    model.setParam('TimeLimit', time_limit)\n",
        "    \n",
        "    print(f\"\\n{'='*80}\")\n",
        "    print(f\"Solving with MIPFocus = {mipfocus_value}\")\n",
        "    print(f\"{'='*80}\\n\")\n",
        "    \n",
        "    # Solve the model\n",
        "    model.optimize()\n",
        "    \n",
        "    # Collect metrics\n",
        "    metrics = {\n",
        "        'MIPFocus': mipfocus_value,\n",
        "        'Status': model.Status,\n",
        "        'StatusName': model.getAttr(GRB.Attr.Status),\n",
        "        'Runtime': model.Runtime,\n",
        "        'NodeCount': model.NodeCount,\n",
        "        'IterCount': model.IterCount,\n",
        "    }\n",
        "    \n",
        "    # Add solution-specific metrics if available\n",
        "    if model.SolCount > 0:\n",
        "        metrics['ObjVal'] = model.ObjVal\n",
        "        metrics['ObjBound'] = model.ObjBound\n",
        "        metrics['MIPGap'] = model.MIPGap\n",
        "        metrics['SolCount'] = model.SolCount\n",
        "    else:\n",
        "        metrics['ObjVal'] = None\n",
        "        metrics['ObjBound'] = model.ObjBound if hasattr(model, 'ObjBound') else None\n",
        "        metrics['MIPGap'] = None\n",
        "        metrics['SolCount'] = 0\n",
        "\n",
        "    \n",
        "    # Clean up\n",
        "    model.dispose()\n",
        "    env.dispose()\n",
        "    \n",
        "    return metrics"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6125dd1c",
      "metadata": {},
      "source": [
        "#### MIPFocus = 0 (Default - Balanced Approach)\n",
        "\n",
        "The default MIPFocus setting strikes a balance between finding new feasible solutions and proving optimality. This is typically where you want to start."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 69,
      "id": "b1d00150",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "================================================================================\n",
            "Solving with MIPFocus = 0\n",
            "================================================================================\n",
            "\n",
            "MIPFocus            : 0\n",
            "Status              : 2\n",
            "Runtime             : 0.2773580551147461\n",
            "NodeCount           : 5816.0\n",
            "IterCount           : 10464.868172922967\n",
            "ObjVal              : 1.0\n",
            "ObjBound            : 1.0\n",
            "MIPGap              : 0.0\n",
            "SolCount            : 10\n"
          ]
        }
      ],
      "source": [
        "# Solve with MIPFocus = 0\n",
        "results_0 = solve_with_mipfocus(mipfocus_value=0)\n",
        "\n",
        "for key, value in results_0.items():\n",
        "    if key != 'StatusName':\n",
        "        print(f\"{key:20s}: {value}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "24133136",
      "metadata": {},
      "source": [
        "#### MIPFocus = 1 (Focus on Feasibility)\n",
        "\n",
        "With MIPFocus = 1, this focuses more on finding good feasible solutions quickly. This is useful when:\n",
        "- You need solutions more quickly\n",
        "- Proving optimality is less critical\n",
        "- The model is difficult and finding any feasible solution is challenging"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 70,
      "id": "18a8c256",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "================================================================================\n",
            "Solving with MIPFocus = 1\n",
            "================================================================================\n",
            "\n",
            "MIPFocus            : 1\n",
            "Status              : 2\n",
            "Runtime             : 0.4498732089996338\n",
            "NodeCount           : 6744.0\n",
            "IterCount           : 11554.081539166284\n",
            "ObjVal              : 1.0\n",
            "ObjBound            : 1.0\n",
            "MIPGap              : 0.0\n",
            "SolCount            : 10\n"
          ]
        }
      ],
      "source": [
        "# Solve with MIPFocus = 1\n",
        "results_1 = solve_with_mipfocus(mipfocus_value=1)\n",
        "\n",
        "for key, value in results_1.items():\n",
        "    if key != 'StatusName':\n",
        "        print(f\"{key:20s}: {value}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1ab7b46d",
      "metadata": {},
      "source": [
        "#### MIPFocus = 2 (Focus on Optimality)\n",
        "\n",
        "With MIPFocus = 2, the focus is more on proving optimality. Try this when:\n",
        "- The solver finds good solutions quickly\n",
        "- The optimality gap closes slowly"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 71,
      "id": "8db07d31",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "================================================================================\n",
            "Solving with MIPFocus = 2\n",
            "================================================================================\n",
            "\n",
            "MIPFocus            : 2\n",
            "Status              : 2\n",
            "Runtime             : 0.23417115211486816\n",
            "NodeCount           : 9766.0\n",
            "IterCount           : 16301.495768467246\n",
            "ObjVal              : 1.0\n",
            "ObjBound            : 1.0\n",
            "MIPGap              : 0.0\n",
            "SolCount            : 10\n"
          ]
        }
      ],
      "source": [
        "# Solve with MIPFocus = 2\n",
        "results_2 = solve_with_mipfocus(mipfocus_value=2)\n",
        "\n",
        "for key, value in results_2.items():\n",
        "    if key != 'StatusName':\n",
        "        print(f\"{key:20s}: {value}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d95a632b",
      "metadata": {},
      "source": [
        "#### MIPFocus = 3 (Focus on Bound)\n",
        "\n",
        "With MIPFocus = 3, we focus on improving the best objective bound. So try when:\n",
        "\n",
        "- The best objective bound is moving very slowly or not at all\n",
        "- You want to make progress on the dual bound (lower bound for minimization, upper bound for maximization)\n",
        "- The solver is struggling to prove optimality because the bound improvement has stalled"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 72,
      "id": "0c802b06",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "================================================================================\n",
            "Solving with MIPFocus = 3\n",
            "================================================================================\n",
            "\n",
            "MIPFocus            : 3\n",
            "Status              : 2\n",
            "Runtime             : 0.19624996185302734\n",
            "NodeCount           : 4358.0\n",
            "IterCount           : 8187.345183761711\n",
            "ObjVal              : 1.0\n",
            "ObjBound            : 1.0\n",
            "MIPGap              : 0.0\n",
            "SolCount            : 10\n"
          ]
        }
      ],
      "source": [
        "# Solve with MIPFocus = 3\n",
        "results_3 = solve_with_mipfocus(mipfocus_value=3)\n",
        "\n",
        "for key, value in results_3.items():\n",
        "    if key != 'StatusName':\n",
        "        print(f\"{key:20s}: {value}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "82569aad",
      "metadata": {},
      "source": [
        "#### Comparative Analysis\n",
        "\n",
        "Now let's compare the results across each of the four MIPFocus values to understand their impact on the solve."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 73,
      "id": "28713f65",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>MIPFocus</th>\n",
              "      <th>Status</th>\n",
              "      <th>Runtime</th>\n",
              "      <th>NodeCount</th>\n",
              "      <th>IterCount</th>\n",
              "      <th>ObjVal</th>\n",
              "      <th>ObjBound</th>\n",
              "      <th>MIPGap</th>\n",
              "      <th>SolCount</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>0</td>\n",
              "      <td>2</td>\n",
              "      <td>0.277358</td>\n",
              "      <td>5816.0</td>\n",
              "      <td>10464.868173</td>\n",
              "      <td>1.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>1</td>\n",
              "      <td>2</td>\n",
              "      <td>0.449873</td>\n",
              "      <td>6744.0</td>\n",
              "      <td>11554.081539</td>\n",
              "      <td>1.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>2</td>\n",
              "      <td>2</td>\n",
              "      <td>0.234171</td>\n",
              "      <td>9766.0</td>\n",
              "      <td>16301.495768</td>\n",
              "      <td>1.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>3</td>\n",
              "      <td>2</td>\n",
              "      <td>0.196250</td>\n",
              "      <td>4358.0</td>\n",
              "      <td>8187.345184</td>\n",
              "      <td>1.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   MIPFocus  Status   Runtime  NodeCount     IterCount  ObjVal  ObjBound  \\\n",
              "0         0       2  0.277358     5816.0  10464.868173     1.0       1.0   \n",
              "1         1       2  0.449873     6744.0  11554.081539     1.0       1.0   \n",
              "2         2       2  0.234171     9766.0  16301.495768     1.0       1.0   \n",
              "3         3       2  0.196250     4358.0   8187.345184     1.0       1.0   \n",
              "\n",
              "   MIPGap  SolCount  \n",
              "0     0.0        10  \n",
              "1     0.0        10  \n",
              "2     0.0        10  \n",
              "3     0.0        10  "
            ]
          },
          "execution_count": 73,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Create a comparison dataframe\n",
        "comparison_df = pd.DataFrame([results_0, results_1, results_2, results_3])\n",
        "comparison_df = comparison_df.drop('StatusName', axis=1)\n",
        "comparison_df"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e7468160",
      "metadata": {},
      "source": [
        "<div style=\"\n",
        "  border-left: 6px solid #3b82f6;  /* blue */\n",
        "  padding:14px 16px;\n",
        "  border-radius:10px;\n",
        "  margin:12px 0;\n",
        "  font-size:16px; line-height:1.4;\">\n",
        "  <strong style=\"font-size:18px;\">🧠 Explore the logs!</strong><br>\n",
        "  Each log file should have saved where you have the mps file. Open them up and start to compare. \n",
        "  \n",
        "  For example, which value of MIPFocus found a better solution faster?\n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "28fa6906",
      "metadata": {},
      "source": [
        "### Quadratic assignment problem"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0f6ced06",
      "metadata": {},
      "source": [
        "Here is a somewhat random variation of the `quadratic assignment problem`. Here is a formulation, but we really won't need it but it's here for fun. "
      ]
    },
    {
      "cell_type": "markdown",
      "id": "accbf59b",
      "metadata": {},
      "source": [
        "#### Mathematical Formulation\n",
        "\n",
        "##### Sets and Indices\n",
        "- $i, k \\in \\{0, 1, \\ldots, n-1\\}$: Facilities\n",
        "- $j, l \\in \\{0, 1, \\ldots, n-1\\}$: Locations\n",
        "\n",
        "##### Parameters\n",
        "- $f_{ik}$: Flow between facility $i$ and facility $k$ (represented as `flow[i][k]`)\n",
        "- $d_{jl}$: Distance between location $j$ and location $l$ (represented as `dist[j][l]`)\n",
        "- $n$: Number of facilities (and locations)\n",
        "\n",
        "##### Decision Variables\n",
        "- $x_{ij} \\in \\{0, 1\\}$: Binary variable equal to 1 if facility $i$ is assigned to location $j$, 0 otherwise\n",
        "\n",
        "##### Objective Function\n",
        "\n",
        "$$\n",
        "\\min \\sum_{i=0}^{n-1} \\sum_{j=0}^{n-1} \\sum_{k=0}^{n-1} \\sum_{l=0}^{n-1} c_{ik} \\cdot f_{ik} \\cdot d_{jl} \\cdot x_{ij} \\cdot x_{kl}\n",
        "$$\n",
        "\n",
        "where the coefficient $c_{ik}$ is defined as:\n",
        "\n",
        "$$\n",
        "c_{ik} = \\begin{cases} \n",
        "-1 & \\text{if } (i+k) \\bmod 2 = 0 \\\\\n",
        "1 & \\text{if } (i+k) \\bmod 2 = 1\n",
        "\\end{cases}\n",
        "$$\n",
        "\n",
        "##### Constraints\n",
        "\n",
        "**Assignment constraints - Each facility to exactly one location:**\n",
        "$$\n",
        "\\sum_{j=0}^{n-1} x_{ij} = 1 \\quad \\forall i \\in \\{0, 1, \\ldots, n-1\\}\n",
        "$$\n",
        "\n",
        "**Assignment constraints - Each location receives exactly one facility:**\n",
        "$$\n",
        "\\sum_{i=0}^{n-1} x_{ij} = 1 \\quad \\forall j \\in \\{0, 1, \\ldots, n-1\\}\n",
        "$$\n",
        "\n",
        "**Binary constraints:**\n",
        "$$\n",
        "x_{ij} \\in \\{0, 1\\} \\quad \\forall i, j \\in \\{0, 1, \\ldots, n-1\\}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9cd14a95",
      "metadata": {},
      "source": [
        "Use this data and the model below to solve the problem. You may notice it take a little while. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 82,
      "id": "766eda9d",
      "metadata": {},
      "outputs": [],
      "source": [
        "n = 10\n",
        "\n",
        "# Distance between locations\n",
        "dist = [[abs(j - k) + 1 for k in range(n)] for j in range(n)]\n",
        "\n",
        "# Flow/weight between facilities\n",
        "flow = [[(i + k + 1) * ((i - k) ** 2 + 1) for k in range(n)] for i in range(n)]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cf80d0f9",
      "metadata": {},
      "source": [
        "Here is the model. You shouldn't have to edit anything now. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 83,
      "id": "a8bc2b06",
      "metadata": {
        "id": "a8bc2b06",
        "outputId": "a0a4b967-14dd-41f3-eabc-f2c1cefff934"
      },
      "outputs": [],
      "source": [
        "# Quadratic assignment problem\n",
        "m1 = gp.Model()\n",
        "\n",
        "# x[i,j] = 1 if facility i is assigned to location j\n",
        "x = m1.addVars(n, n, vtype=GRB.BINARY, name=\"x\")\n",
        "\n",
        "# Build objective function\n",
        "obj = gp.QuadExpr()\n",
        "for i in range(n):\n",
        "    for j in range(n):\n",
        "        for k in range(n):\n",
        "            for l in range(n):\n",
        "                coeff = flow[i][k] * dist[j][l]\n",
        "                # Alternate signs\n",
        "                if (i + k) % 2 == 0:\n",
        "                    coeff = -coeff\n",
        "                obj.addTerms(coeff, x[i, j], x[k, l])\n",
        "\n",
        "m1.setObjective(obj, GRB.MINIMIZE)\n",
        "\n",
        "# Each facility is assigned to exactly one location\n",
        "for i in range(n):\n",
        "    m1.addConstr(x.sum(i, \"*\") == 1, name=f\"assign_facility[{i}]\")\n",
        "\n",
        "# Each location is assigned exactly one facility\n",
        "for j in range(n):\n",
        "    m1.addConstr(x.sum(\"*\", j) == 1, name=f\"assign_to_location[{j}]\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ef36ffdc",
      "metadata": {},
      "source": [
        "Run this without editing any parameters to see how long it takes. Then, try to recall a parameter mentioned in the log file session that may fit this type of problem. Try searching the [documentation](https://docs.gurobi.com/projects/optimizer/en/current/).\n",
        "\n",
        "- Hint 1: Look at the type of terms you see in the log file and think about how these terms are typically dealt with. \n",
        "- Hint 2: We tend to like things *linear*.\n",
        "- Hint 3: You'll be told what it is in a couple of minutes. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 84,
      "id": "1b7940cf",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[rosetta2] - Darwin 24.6.0 24G325)\n",
            "\n",
            "CPU model: Apple M4\n",
            "Thread count: 10 physical cores, 10 logical processors, using up to 10 threads\n",
            "\n",
            "Optimize a model with 20 rows, 100 columns and 200 nonzeros (Min)\n",
            "Model fingerprint: 0xb35f04ea\n",
            "Model has 0 linear objective coefficients\n",
            "Model has 5050 quadratic objective terms\n",
            "Variable types: 0 continuous, 100 integer (100 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 1e+00]\n",
            "  Objective range  [0e+00, 0e+00]\n",
            "  QObjective range [2e+00, 3e+04]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [1e+00, 1e+00]\n",
            "Found heuristic solution: objective 1644.0000000\n",
            "Presolve time: 0.02s\n",
            "Presolved: 20 rows, 100 columns, 200 nonzeros\n",
            "Presolved model has 5050 quadratic objective terms\n",
            "Variable types: 0 continuous, 100 integer (100 binary)\n",
            "\n",
            "Root relaxation: objective -6.799610e+05, 119 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 -420272.91    0  100 1644.00000 -420272.91      -     -    0s\n",
            "H    0     0                    -16804.00000 -420272.91  2401%     -    0s\n",
            "H    0     0                    -18856.00000 -420272.91  2129%     -    0s\n",
            "H    0     0                    -19144.00000 -420272.91  2095%     -    0s\n",
            "H    0     0                    -25848.00000 -420272.91  1526%     -    0s\n",
            "H    0     0                    -25972.00000 -420272.91  1518%     -    0s\n",
            "H    0     0                    -26036.00000 -420272.91  1514%     -    0s\n",
            "     0     0 -420272.91    0  100 -26036.000 -420272.91  1514%     -    0s\n",
            "H    0     0                    -28500.00000 -420272.91  1375%     -    0s\n",
            "H    0     0                    -28788.00000 -420272.91  1360%     -    0s\n",
            "     0     2 -420272.91    0  100 -28788.000 -420272.91  1360%     -    0s\n",
            "H  190   224                    -29508.00000 -420272.91  1324%   3.5    0s\n",
            " 286876 77958 -32763.711   44    4 -29508.000 -186216.59   531%   3.4    5s\n",
            " 629948 111790 -40560.062   45    4 -29508.000 -141523.49   380%   3.4   10s\n",
            " 984187 114874 -112738.58   38   13 -29508.000 -125830.24   326%   3.4   15s\n",
            " 1342100 97214 -52887.370   47    4 -29508.000 -98812.913   235%   3.4   20s\n",
            " 1688047 21840 -29644.586   34    4 -29508.000 -67628.636   129%   3.4   25s\n",
            "\n",
            "Explored 1722990 nodes (5879798 simplex iterations) in 25.54 seconds (38.70 work units)\n",
            "Thread count was 10 (of 10 available processors)\n",
            "\n",
            "Solution count 10: -29508 -28788 -28500 ... 1644\n",
            "No other solutions better than -29508\n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective -2.950800000000e+04, best bound -2.950800000000e+04, gap 0.0000%\n"
          ]
        }
      ],
      "source": [
        "# Set parameters here\n",
        "# m.Params.ParameterName = value\n",
        "m1.optimize()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "108f1733",
      "metadata": {},
      "source": [
        "Here's the whole model with the updated parameter. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 77,
      "id": "151cfb80",
      "metadata": {
        "id": "151cfb80",
        "outputId": "9dc1798b-6050-4885-8db2-7f8830c0769b"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter PreQLinearize to value 2\n",
            "Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (mac64[rosetta2] - Darwin 24.6.0 24G325)\n",
            "\n",
            "CPU model: Apple M4\n",
            "Thread count: 10 physical cores, 10 logical processors, using up to 10 threads\n",
            "\n",
            "Non-default parameters:\n",
            "PreQLinearize  2\n",
            "\n",
            "Optimize a model with 18 rows, 81 columns and 162 nonzeros (Min)\n",
            "Model fingerprint: 0x05cc9672\n",
            "Model has 0 linear objective coefficients\n",
            "Model has 3321 quadratic objective terms\n",
            "Variable types: 0 continuous, 81 integer (81 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 1e+00]\n",
            "  Objective range  [0e+00, 0e+00]\n",
            "  QObjective range [2e+00, 2e+04]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [1e+00, 1e+00]\n",
            "Found heuristic solution: objective 765.0000000\n",
            "Presolve time: 0.00s\n",
            "Presolved: 98 rows, 161 columns, 3562 nonzeros\n",
            "Variable types: 0 continuous, 161 integer (82 binary)\n",
            "\n",
            "Root relaxation: objective -2.007734e+05, 116 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 -200773.38    0   64  765.00000 -200773.38      -     -    0s\n",
            "H    0     0                    -945.0000000 -200773.38      -     -    0s\n",
            "H    0     0                    -2701.000000 -200773.38  7333%     -    0s\n",
            "H    0     0                    -8903.000000 -163034.43  1731%     -    0s\n",
            "H    0     0                    -9401.000000 -163034.43  1634%     -    0s\n",
            "     0     0 -163034.43    0   94 -9401.0000 -163034.43  1634%     -    0s\n",
            "H    0     0                    -13073.00000 -158309.29  1111%     -    0s\n",
            "     0     0 -158309.29    0  103 -13073.000 -158309.29  1111%     -    0s\n",
            "     0     0 -158052.61    0  104 -13073.000 -158052.61  1109%     -    0s\n",
            "     0     0 -131037.07    0  118 -13073.000 -131037.07   902%     -    0s\n",
            "     0     0 -130602.04    0  118 -13073.000 -130602.04   899%     -    0s\n",
            "     0     2 -130602.04    0  118 -13073.000 -130602.04   899%     -    0s\n",
            "H   80   122                    -15377.00000 -130602.04   749%  26.1    0s\n",
            "H   85   122                    -16721.00000 -130602.04   681%  26.5    0s\n",
            "H  329   272                    -18433.00000 -113575.06   516%  21.3    0s\n",
            "H  423   371                    -19025.00000 -113575.06   497%  20.7    0s\n",
            "H 1160   669                    -19219.00000 -83118.181   332%  23.9    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Gomory: 1\n",
            "  Implied bound: 78\n",
            "  MIR: 155\n",
            "  StrongCG: 2\n",
            "  Inf proof: 2\n",
            "  Zero half: 1\n",
            "  RLT: 61\n",
            "\n",
            "Explored 9726 nodes (209583 simplex iterations) in 0.93 seconds (1.41 work units)\n",
            "Thread count was 10 (of 10 available processors)\n",
            "\n",
            "Solution count 10: -19219 -19025 -18433 ... -945\n",
            "No other solutions better than -19219\n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective -1.921900000000e+04, best bound -1.921900000000e+04, gap 0.0000%\n"
          ]
        }
      ],
      "source": [
        "n = 9\n",
        "\n",
        "# Quadratic assignment problem\n",
        "m2 = gp.Model()\n",
        "# Distance between locations\n",
        "dist = [[abs(j - k) + 1 for k in range(n)] for j in range(n)]\n",
        "\n",
        "# Flow/weight between facilities\n",
        "flow = [[(i + k + 1) * ((i - k) ** 2 + 1) for k in range(n)] for i in range(n)]\n",
        "\n",
        "# x[i,j] = 1 if facility i is assigned to location j\n",
        "x = m2.addVars(n, n, vtype=GRB.BINARY, name=\"x\")\n",
        "\n",
        "# Build objective function\n",
        "obj = gp.QuadExpr()\n",
        "for i in range(n):\n",
        "    for j in range(n):\n",
        "        for k in range(n):\n",
        "            for l in range(n):\n",
        "                coeff = flow[i][k] * dist[j][l]\n",
        "                # Alternate signs\n",
        "                if (i + k) % 2 == 0:\n",
        "                    coeff = -coeff\n",
        "                obj.addTerms(coeff, x[i, j], x[k, l])\n",
        "\n",
        "m2.setObjective(obj, GRB.MINIMIZE)\n",
        "\n",
        "# Each facility is assigned to exactly one location\n",
        "for i in range(n):\n",
        "    m2.addConstr(x.sum(i, \"*\") == 1, name=f\"assign_facility[{i}]\")\n",
        "\n",
        "# Each location is assigned exactly one facility\n",
        "for j in range(n):\n",
        "    m2.addConstr(x.sum(\"*\", j) == 1, name=f\"assign_to_location[{j}]\")\n",
        "\n",
        "# Set parameters here\n",
        "m2.Params.PreQLinearize = 2\n",
        "\n",
        "m2.optimize()"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "gurobi_ml",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.11.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}
